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Abstract 

We analyze the biquadratic bilinear Heisenberg magnet on a hon- 
Q_f eycomb lattice via Schwinger boson formalism. Due to their vulner- 

£^ \ ability to quantum fluctuations, non conventional lattices (kagome, 

triangular and honeycomb for example) have been cited as candidates 
to support spin liquid states. Such states without long range order 
at zero temperature are known in one-dimensional spin models but 
their existence in higher dimensional systems is still under debate. Bi- 
quadratic interaction is responsible to various possibilities and phases 
as it is well-founded for one-dimensional systems. Here, we have used a 
bosonic representation to study the properties at zero and at finite low 
temperatures of the biquadratic term in the two-dimensional hexag- 
onal honeycomb lattice. The results shown a ordered state at zero 
temperature but much more fragile than that of a square lattice; the 
behavior at finite low temperatures is in accordance with expectations. 

1 Introduction 

Non conventional lattices in magnetic systems have been received much 
attention in recent years. Traditional square lattices are well established and 
no surprises are expected. On the other hand, the non-conventional lattices 
are serious candidates to the so-called two-dimensional spin liquid phase. 
Spin liquids are disordered states of matter with power-law decay of spin- 
spin correlations and zero local magnetic moment. That states occur at zero 
temperature and their disorder is derived from quantum and not thermal 
fluctuations. Such phases are known to exist in one- dimensional antiferro- 
magnets (AF) but, in higher dimensions, they still remain uncertain. Most 
known two-dimensional magnetic materials have a Neel order at zero temper- 
ature (T = 0) even though, some unusual systems are candidates to present 
a spin liquid state. In three-dimensional magnetic systems, the existence of 
this state is even more unexpected due to stronger spin interactions. 

One possible way to obtain the 2D spin liquid is through geometric frus- 
tration presents in some lattices. In the classical antiferromagnetic kagome 
lattice, for example, it is impossible to align all neighbors spins and the 
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ground state is highly degenerated. Anisotropics together with longer range 
interactions (second and far nearest neighbors) contribute to disorder the 
ground state and to increase the possibility of a spin liquid state. The prop- 
erties at zero temperature are questionable but recent works indicate the 
occurrence of a spin liquid state [H|3j|lJ|6] • Other lattices, although non frus- 
trated, have been shown interesting possibilities as it is the case of hexagonal 
honeycomb lattice. The honeycomb is the two-dimensional case with the 
smallest coordination number (neighbors number) z — 3. It is between the 
disordered one dimensional spin model with z = 2 and the ordered (in zero 
temperature) square lattice with z — 4. Thus such a system may have 
larger vulnerability to quantum fluctuations, mainly for small spins, and it 
is a possible candidate to be a two-dimensional spin liquid . Nevertheless, 
recent works have shown a Neel order for the spin- 1/2 Heisenberg AF on 
honeycomb lattice j5j|6], although weaker than the square lattice case; In 
addition, frustrated honeycomb models have been demonstrated disordered 
ground states HHIIJD1]. 

In present work we are interested in the behavior of the biquadratic bi- 
linear Heisenberg model in a honeycomb hexagonal lattice. It is well known 
that the Heisenberg Hamiltonian describes both the ferromagnetic and anti- 
ferromagnetic materials in any dimension. Here we add the biquadratic term 
in such a way that we study the following model: 



where the sum is over nearest neighbors and constants J\ and J2 define the 
bilinear and biquadratic interactions, respectively. It is important to high- 
light that we adopting 5=1 once the biquadratic term makes sense only to 
S > 1/2. The SU(2) spin rotation symmetry in Hamiltonian ([T]) is preserved 
and so we can expect Goldstone modes as lower energy excitations over the 
ground state. We consider the case Ji = 1 (antiferromagnetic coupling) and 
— J\ < J2 < Ji, but it is common to write J\ = cos# and J2 = sind. The 
one-dimensional case is well documented [T2TT14"] and the various phases were 
already discovered. For 9 = tt one has the usual ferromagnet while 9 = 
corresponds to the usual Heisenberg antiferromagnet system. On the range 
71 < 9 < 5tt/4, there is a stable ferromagnetic regime with long range order 
(LRO); for — tt/A < 9 < 7r/4 there is an antiferromagnetic phase with Hal- 
dane gap (spin-1) and, in the limit n/49 < tt/2, there is a trimerized phase. 
Some points have exact solution. For instance, the angles 9 = ±7r/4 can be 
solved by the Bethe ansatz and beyond the angle tan# = 1/3 corresponds to 
AKLT Model. The two-dimensional case is more complicated and only some 
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regions are known. Ivanov et al. [T^jlTE] have shown a nematic phase for 
9 > 5ir/4 by using a continuous model similar to nonlinear sigma model (the 
same result has been achieved by Chubukov using the Holstein-Primakoff 
representation [H|); for 5tt/5 < 9 < 77r/4 there is a dimerized ferromagnetic 
phase. 

We can represent the biquadratic term as a function of spin and quadrupole 
operators (second-order spin moment): 



,2 1 ,~ - , 1 ,~ - s 4 



(S« • SjY = - (Qi ■ Qj) --{St- Sj) + -, (2) 
where Qi operators are given by: 

g <P) = 2(gf) 2 -(gf) 2 -(gy) 2 > (3) 

V3 

Qf ) = (Sf) 2 - (ST) 2 , (4) 

= Sf£% + SfSf. (5) 
Qf = 5?5? + 5?5?, (6) 
Qf = 5?5f + 5?5?. (7) 

The three spin operators together the five quadrupole operators form the 
generators of the SU(3) Lie group. For the special case Jl = J2 one has 
the SU(3) symmetric ferromagnetic model [IB] while for Jl = one has the 
SU(3) symmetric valence-bound antiferromagnet |l9fl2Tj . 

The methods used to study systems in condensed matter physics, mag- 
netism especially, are vast and diversified. In the current work, we have 
used the Schwinger boson representation to develop the physics of the bi- 
linear biquadratic Heisenberg model in a honeycomb lattice at both, zero e 
low temperatures. The Schwinger formalism has some advantages over other 
bosonic representations (Holstein-Primakoff and Dyson-Maleev). The first is 
the holonomic constraint that fix the number of bosons and can be imple- 
mented easily by a Lagrange multiplier. Second, there are not root terms 
and so we do not need to specify a preferential direction to the ground state, 
as occur in the Holstein-Primakoff method. Therefore we are able to treat 
ordered and disordered phases what are important in search for a spin liquid 
state. Following the usual prospects we have adopted the boson condensa- 
tion at zero temperature to avoid the divergences of the theory and, at low 
temperatures, we have used approximations that, within the correct limits, 
provide coherent results. The paper is as following: in Sec. |2] we devel- 
oped the Schwinger bosons formalism; in Sec. E] we present the results for 
zero temperature and finite low temperature and, finally, the conclusions are 
exposed in last section (j2j). 
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2 Formalism 



Commonly, the spin operators are defined by two 577(2) Schwinger opera- 
tors but, because the biquadratic term, we have considered the representation 
by S77(3) Schwinger formalism. Thus, each spin operator is represented by 
three bosonic operators a irn , where i denotes the site % and m = — f, 0, 1 
specifies the eigenvalues of Sf. Accordingly, a| m |0) creates a particle with 
^-component of spin S z = m on site % (|0) is the vacuum in the Fock space). 
The generators F™ n = a\ m a itn form the SU (3) Lie algebra and obey the com- 
mutation relation [F t mn , Ff] = 5ij(5 n:P F™ q - 5 miP F? q ). As a function of the 

a operators, the spin operators are expressed by Sf = >/2(a|o a i,-i + a i l a,i ! o), 
S~ = y/2(a\ idifi + a| aj ; i) and S? = aj^a^i — a\ ^i-i. The bosonic oper- 
ators keep the spin commutation relation and to fix Sf = S(S + 1) we have 
to impose the local constraint X] m a !m a *,m = m or der to symmetrize 
the spin and quadrupole operators, we apply a rotation over the a operators 
defining new operators b as follow: 

hi = ^= (di ,-i - a»,i) , (8) 

bi2 = ^= (<ii ,-i + a it i) , (9) 
b i3 = a ifi . (10) 
Therefore the spin operators are written as: 

= -i (b\ 2 b i3 - b\ z b i2 ) , (11) 

S! = ~i (fei - ftJi&is) , (12) 

SI = ~i (blb i2 - b\ 2 b a ) , (13) 

and the quadrupoles: 

Q<°> = -L (b\ lbil + b \ 2 b i2 - 2b\ 3 b l3 ) , (14) 

Qt ) = - (frJibii - 4 2 b i2 ) , (15) 

<r = - (fe + , (is) 

QT = - (b\ 2 h 3 + b\,b i2 ) , (17) 
QT = ~ (fei + blhs) , (18) 
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while the constraint holds the same. The biquadritic bilinear Heisenberg 
hamiltonian as a function of b operators is expressed by: 

H = E - J ±)AAi + J i : 4% : ] > ( 19 ) 

(<j> 

where we have introduced the bond operators A%j= ^„ ^i/j^jn an d &ij = 
S«%A>> ^ ne ^ wo P om ts denote normal ordered. The hamiltonian is fourth 
order in b and we decouple it by using the Hubbard- Stratonovich [221123] 
transform: 

Ej^j (4>Sy + (E^Hj. - (Sj,)<Hy>. (20) 

In above equation Sy = Aij,Bij, where we have adopted A = (A\j) = (Aij) 
and B = {B\-) = (Bij). Therefore the second order mean field hamiltonian 
is given by: 

H MF = [( J 2 - J,) A 2 + JiB 2 } -NSX + XJ2 



2 

[(•& - -AM (4 + A,) + JiS (4 + fit 



(21^ 



The constraint ^l^A/t — 5 1 was implemented by a Lagrange multiplier Aj 
on each site and we have also adopted a constant value A = (Aj). The mean 
field parameters A, B and A are determined by minimizing the Helmholtz 
free energy. 

The honeycomb is a bipartite but not a Bravais lattice and so we have 
to treat each sublattice separately. The sublattices R and R' are hexagonal 
Bravais lattices and they are coupled by nearest neighbors interactions. After 
Fourier transforming the Schwinger bosons independently on each sublattice, 
we obtain 

k 

and: 



N 

k 



the Hamiltonian is written as: 



k /i 

+3 (^/ k ,7k + b'^l) + A (bib + JJ.C. , (24) 
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where H are constant terms and 7k = e l¥>k |7k| is the structure factor: 



7 k =3 



\f3ky 



cos — cos 
2 2 



+ cosk x + 2isen-^-cos 



y/3k y 



isenk* 



(25) 



In the Hamiltonian (12~3|) . the two sublattices are still coupled and dif- 
ferently of the square lattice, the structure factor for honeycomb lattice 
is not real. We resolve these two problems defining new operators 6 kA1 = 

^jr ( ic L„ + 41) and = (~ ic ^ + c ki)- The new bosons c k hold 

all commutation relations and the Hamiltonian (as a function of them) is 
written as: 



k 



n 

k • 



(26) 



with f3^ = (c kl , c k 2, c k tj, ci kl , cf_ k2 , ci k3 ) where s = /, If, while the matrices 
are H 1 = (\ - 3JiB|Tk|)<7„ ® 4x3 + 3(J 2 - Ji)A|7 k |a x <g> J 3X 3 and = 
(A + 3Ji-B|7 k |)cr <g> I 3 x3 + 3(J 2 - ^i)>l|7k|^ <8> 4x3 (here <7; are the Pauli 
matrices). if MF can be diagonalized by a canonical Bogoliubov transform: 



c kjU = cosh^a^ + senh^c/"^ 
c£ = cosh^ag, + senh# k V 



(27) 
(28) 



We choose and dy. 1 so that the non-diagonal terms vanish, what is done 
by: 



tanh2#£ 



tanh2^ 



A-3J!ff|7 k | 
'3(J 2 -^i)A|7k| 
A + 3J!5|7 k | 



3(J 2 - Ji)A| 7k | ' 
Once diagonalized, H MF gives the eigenvalues of energy: 



(29) 
(30) 



(A - ZJxB^f - (3A(J 2 - Ji)| 7 k|) 2 
(A + 3J lJ B| 7k |) 2 -(3A(J 2 - Ja)| 7 k|) 2 . 



(31) 
(32) 



To determine the parameters A, B and A we need to minimize the Helmholtz 
free energy: 



senh 



+ In 



senh 



0E n 



(33) 



what yields the mean field self-consistent equations: 
3 3 



A = 

and 
B = 



2N 



E 



(()tli f3E^ A-3J!g| 7 k| | cQth f/3E n \ A + 3JiS|7 k | 



E 
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2N 



3 

2iV 



A + 3Ji5|7 k | 



(34) 

V [coth m\ 3 ^" Jl) + coth f ^ 3 ^" Jl) l l7 k | 2 

k L V 2 / E i V 2 / % J 

(35) 

, — ^ I- — i-i,ki //SE/N A-3JiB|7k|] . . 
coth I — — I coth I — ^ — 1 — |7k|. 

(36) 

As well known, one and two-dimensional systems can have LRO only at zero 
temperature and this implies an abrupt chance at T = 0. Indeed when 
the temperature approaches the zero, the bosons condensate at a zero en- 
ergy state and the self-consistent equations diverge. This inconvenience is 
surmounted by separating the divergent term of the sum and introducing a 
new parameter (the condensate density) as it is done in the Bose-Einstein 
condensate. For finite temperatures that problem does not exist and the 
equations can be resolved normally. Obviously the self-consistent equation 
can not be solved exactly and numeric methods or approximations have to 
be used. In the next section we present the results for both T = and low 
finite temperatures. 

Using the equations (|2"7|) and (12 8| . we calculate the mean field double 
boson operators: 



1 



k 



V 1 g-ikAr 



cosh20£ 



1 



cosh20£ J 



ji 



ikAr 



-sinh20£ ( n{ + J ) + sinh2^ J [r£ + - 



(37) 
(38) 



for i and j belonging to the same sublattice and: 



-cosh2^ ( nl + = ) + cosh2C ( nl 1 



sinh2^ ( n{ + J ) + sinh2^ 7 ( n£ J + - 



(39) 



(40) 
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for i and j of different sublattices. The phase angle (p^ is the same that 
appears in the structure factor (equation (1251) ). whilst the boson densities 
are given by: 

< = («) = ( 41 ) 

with t — I, II. All others mean field are null. 

3 Results 

When the temperature decreases to zero, a transition phase takes a place 
and one of the energies of the spectrum vanishes, characterizing a boson con- 
densation. Therefore the self-consistent equations diverge and there are not 
more solutions of the parameters A, B and A. According to Takahashi and 
Arovas et al. |24T - [2~T] . the non-existence of solutions is related to a sponta- 
neous broken symmetry, since Schwinger formalism is invariant over SU(2). 
At finite temperatures, there are solutions for any dimension and the system 
is disordered, i.e., there is not long range order (at zero temperature, there are 
solutions only for the one-dimensional case). Here, the condensation occurs 
for Ej if Ji > and En otherwise. The ground state is therefore ordered and 
the lowest excitation energies are the massless Goldstone modes. Expanding 
for small k, the dispersion relations (|3Tj) and fl32|) assume a relativist form: 



E I = ^A] + k 2 c 2 and E n = ^ A 2 H + k 2 c 2 /; (42) 

with the gap energies: 

A/ = y/(X - ^hBj - (3( Ji - J 2 )A)\ (43) 

A n = y/(X + 3JiB)2 - (3(Jx - J 2 )A)2, (44) 

while the spin-wave velocities are given by: 

ci = sj\ (3JiB - 9J*B* + 9(Ji - J 2 )M 2 ) (45) 

, c H = (-37x5 - 9J 2 5 2 + 9(Jx - J 2 )M2). (46) 



The spin-wave velocities as a function of J 2 (Ji = 1) are plotted in figure [TJ 
A similar behavior is observed for the frustrated honeycomb Heisenberg [8], 
where there is a linear decreasing of the spin-wave velocity as a function 
of the second-nearest neighbors exchange coupling. For J 2 > 0.65, the c 2 
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spin-wave velocity is null whilst the point where C\ = is beyond of the limits 
considered. 

Considering positive values of J\, the condensation occurs for Ej and 
then, Ai = while An is finite. After separating the divergent terms, the 
self-consistent equations are written in the continuous limit as: 



s 4 



3 f d 2 k IX - 3Jii?|7 k | A + 3JiB\j k \ 
2 / ~2a ~E, + E 



n 



A = p- 



d 2 k 

2a 



3A(J 2 - JO^kJ 2 3A(J 2 - J0| 7k 



21 



E 



ii 



and 



B = -p + 



d 2 k r (A-3J lJ B| 7k |)^ k ] + (A_+3J lJ B| 7k |)| 7k | 
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(47) 
(48) 

(49) 



where p is a new parameter that measures the boson condensate and a = 
is the first Brillouin zone area. 




0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 



-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 



Figure 1: The spin- wave velocities (left) and the boson density condensate 
(right) in function of J 2 



The condensate density is plotted in figure [TJ For J 2 = we have the 
pure Heisenberg antiferromagnet, in which the condensate density is approx- 
imately 0.64. As expected p is smaller than the corresponding one for the 
square lattice (which is approximately 80% smaller [27J). The condensate 
density increases as of J 2 increases (with an almost linear behavior in the 
range —J\ < J 2 < 0.2Ji). Comparing the two graphics in Fig.([T]) we can 
see that the increasing condensate density (and consequently the ordering) 
occurs together with the decreasing of the spin-wave velocities. This is ex- 
pected since the spin-wave is responsible for disordering the ground state and 
then, the higher velocity implies in higher mess. Therefore the honeycomb 
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lattice holds a long range order at T = but due the small coordination 
number this ordering is weaker than that of the square case. 

As it is well known, the two-dimensional square lattice presents LRO 
for all spin values S > S c ~ 0.19. Writing the magnetization as m s = 

S + I - I I ( A ~ 3 gf l7kl + A+3 ^f l7kl ) we determine the critic spin S c for 
which m s — > 0. The results are shown in Fig. [2j The S c curve separates 
the region with an ordered ground state from the disordered ground state. 
For the square case the disordered ground state is inaccessible to all physical 
spins. Using equations ( I37p and l3"Sj) we obtain the mean value (Sj • Sj) ~ 2.5 



Ordered ground state 



Disordered ground state 



-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 



Figure 2: The critic spin value separates the region with a ordered ground 
state from that with a disordered ground state. 

for all values of J2- It is greater than the expected value S(S + 1) = 2 by 
a factor of approximately 3/2. Such factor also appears in the equations 
obtained by AA [25] and this discrepancy arises because we imposed the 
constraint only on the average. The Fourier transform S(q) of the two-point 
function (Sj • Sj) is as follow: 



S(q) = -\ + | E { cosh(2^ +q - 261) (< +q + \) (n£ + \ 



+ 



Hi nam I n n_ + t 



n 



ii 



+ coah(20ft q -20£, r , k+q 2 
+ [1 _ c «w«,-¥*] cosh(2^ +q + 2^) (< +q + ~ 
+ cosh(20& q + 20£) K J +q +^| I//, 



[1 + 



>Y>k+q-¥>k 



] + 



1 



^ + 2> 



(50) 
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The static uniform susceptibility \ 

2 



5(0) 
3T 



is therefore: 



X 



T 



N 



nl+l) 



+ nl I (nl I + l 



(51) 

] for an antiferromagnetic 



Similar equations were obtained by Takahashi 
system in a square lattice. 

For T^O the energies (13 ip and (132 p are finite and we have no divergences 
in the self-consistent equations. The equations (J33J, (I3"5"j) and (I55|) . therefore, 
can be solved using numeric methods but we have adopted an approximation 
following Yoshida [28J. Matching equations (I34p and ( H7|) we got: 

d 2 k 



2a 
d 2 k 
~2V 



A - 3JiB 7k 



A - 3Ji5 |7k| 



coth 



A + 3J 1 5| 7k | 



coth 



/3£ 



A - 3Ji5 |7k| 



, , • ( 52 ) 

where the "0" index indicates the solutions at zero temperature. The above 
equation is solved at the limit of low temperatures. We separate the integral 
in two regions: the first one is a circle around the origin of radius ku and 
the other is the remaining area of the Brioullin zone. The radius is chosen 
such that T <C cku (&b = 1) with the spin-wave velocity c of the same 
order of c\ and c<i- In principle, the superior value for ku is much smaller 
than 1 and it should be chosen so that the energies can be approached by 
Ej = y/tfj + k 2 c 2 and E n = ^/A 2 n + k 2 c 2 H . For fcjvf = 1 the error between 
the exact and approximated energy is around 10%, which allow us to assign 
ku ~ 0.1 as a reasonable limit. Thus, we can estimate ku for not so small 
temperatures. At zero temperature, we have E j = kcoj (massless gap 

mode) and E q ji = \J 11 + ^ 2c on- Inside the first region ! ^| LL w 1 

and so coth (^f -) ~ 1 whilst coth (^§ L ) holds without more approximations. 
In the extern region, k > kuj the energies are not so small to be approached 
by the relativistic dispersion relation; meanwhile the cotangent terms (at 
low temperatures) are taken as unitary. With these considerations and after 
some work, the density condensate is calculated as: 



P 



6tt(A - 3JiS)T 



ac, 



In senh 



ygj + k 2 /C 2 
2T 



In senh ( — ^ ) 
\2T) 



3tt(A - 3Ji J B )k 



M 



(53) 



which, at low temperature limit, yields A/(T) = Te kT with k — 6 _. x _., , ,,, . 
In Fig. ([3|) we show the gap energy Aj as a function of the temperature 
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for three different exchange constants J 2 and hi Aj as a function of J 2 for 
T = 0.15. The finite gap result agrees with the Mermin- Wagner theorem. 




Figure 3: The gap energy Aj as function of temperature for some values of 
J 2 (left) and In Aj as function of J 2 for T = 0.15 (right). 



4 Conclusions 

Using the Schwinger boson formalism we have studied the bilinear bi- 
quadratic Heisenberg model at zero and low temperatures. We have shown 
that, inside the considered limits of J 2 (the biquadratic coupling), the ground 
state at zero temperature remains ordered. Considering the boson conden- 
sation, we have shown that the "degree" of order in a honeycomb lattice is 
between 58% (for J 2 = — 1) and 78% (for J 2 = 1), which is weaker than 
the 81% observed for the square lattice (as expected). Therefore, even with 
a smaller coordination number (z = 3), the quantum fluctuations in the 
honeycomb lattice are not sufficient to create a two-dimensional spin liquid 
state. Our approach is not appropriate for | J 2 | > 1; however, the asymptotic 
behavior of the negative values of J 2 (Figure ([1]) suggests the absence of a dis- 
ordered phase in the limit J 2 — 1 while for J 2 > 1 the system seems to be 
strogly ordered. We have also shown that the ordered ground state exists for 
all physical spin values. The superior value of spin to a phase with (m) = 
is around 0.42 and it occurs when J 2 = —1. Above the ground state, the low 
energy excitations are massless Goldstone modes with relativist dispersion 
relation, since there is a spontaneous broken symmetry. The spin-wave veloc- 
ities decrease almost linearly as a function of J 2 (Figure ([T]) and c 2 vanishes 
when J 2 rs 0.65 ( C\ vanishes for J 2 > 1, outside the limit considered). An- 
alyzing the condensate density graphic, one can see that the slow spin- wave 
corresponds to a more ordered system (higher condesantion). Although the 
Schwinger formalism is not the best way for treating finite temperatures, we 
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have found consistent results. At low temperatures, the ground state is dis- 
ordered and the excitations have a gap that increases with the temperature 
(A oc Te~ KT ) as dictated by the Mermin- Wagner theorem. 
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